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Abstract 

Cortical spreading depression (CSD) is a slow-moving ionic and metabolic disturbance that propagates in 
cortical brain tissue. In addition to massive cellular depolarization, CSD also involves significant changes 
in perfusion and metabolism - aspects of CSD that had not been modeled and are important clinically. In 
this study, we develop a mathematical model for CSD where the sodium-potassium-ATPase, responsible 
for neuronal polarization and recovery from CSD, operates at a rate that is dependent on the supply 
of oxygen. The supply of oxygen is determined by modeling blood flow through a lumped vascular tree 
with an effective local vessel radius that is controlled by the extracellular potassium concentration. Our 
model replicates the qualitative and quantitative behavior of CSD found in experimental studies and 
elucidates the effect of oxygen deprivation on CSD recovery. Wc show that during CSD, the metabolic 
activity of the cortex outstrips the physiological limits placed on oxygen delivery, regardless of vascular 
constriction or dilation. However, vasoconstriction and vasodilation play important roles in determining 
the susceptibility of the cortical tissue to CSD and its recovery. Our study incorporates relevant perfusion 
and metabolic factors into a model of CSD, and in doing so, helps explain the varied phenotypes of this 
phenomenon in vitro and in vivo. 



Author Summary 

Cortical Spreading Depression (CSD) is a traveling wave of electrical disturbance in the brain that is 
associated with migraine and brain injury. Massive changes in cerebral blood flow have been associated 
with CSD, but thus far no mathematical modeling has incorporated these physiologically important 
changes. In this manuscript, wc present a mathematical model for CSD that incorporates blood flow and 
oxygen consumption in order to examine the effects of the brain vasculature on CSD and vice versa. 
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1 Introduction 

Cortical spreading depression (CSD) is a propagated depolarization that occurs in the gray matter of many 
species PQ. In humans, it is known to occur during brain injury, stroke, and subarachnoid hemorrhage [2|. 
There is also strong evidence that CSD is responsible for migraine aura [3H5]. 

Although CSD was discovered in 1944 by Leao [6], we still do not have a detailed understanding of 
how CSD is manifest. In particular, CSD has been associated with massive changes in cortical perfusion. 
The magnitudes of these changes vary by animal species, but significant decreases and increases in blood 
flow or volume occur in all species tested [7HS]. Also common to all species tested is a mismatch in the 
delivery of substrates to meet metabolic demands, i.e., a derangement of neurovascular coupling |1Q[|11) . 

Blood delivery is known to play a significant role in CSD in several ways. Changes in perfusion can 
induce CSD during stroke, called a peri-infarct depolarizations (PID). Conditions that mimic the effects 
of hypoperfusion, such as oxygen glucose deprivation (OGD) and exposure to ouabain (an inhibitor of 
the sodium potassium ATPase), also generate spreading depolarizations [T2HT4"] . 

Perfusion changes can also modulate the signature characteristics of CSD. Depending on the levels 
of the underlying oxygenation or blood pressure, the amplitude and duration of depolarization and the 
velocity of propagation of CSD can be altered [T4TT5] . Clearly, CSD in vivo cannot be understood without 
reference to the vascular changes that condition - and are conditioned by - the phenomenon. There is a 
need to further explore the implications of the effects of perfusion and metabolism on various aspects of 
CSD. 

Such explorations naturally lead to the development of mathematical models in which many mecha- 
nisms can be studied independently and/or simultaneously. Mathematical models of CSD have accounted 
for ionic diffusion, cellular membrane ionic currents, sodium-potassium-ATPasc, and extra- and intra- 
cellular volume changes. Thus far, mathematical models of CSD have not looked at the dynamical 
implications of neurovascular coupling and metabolism. Here we develop a mathematical model that 
uses known physiological data relating effective blood vessel diameter and extracellular potassium con- 
centration to model oxygen delivery in the brain. We show that the oxygen deprivation that results from 
both metabolic demand and vasoconstriction modifies the characteristics of CSD waves and helps explain 
some of the variability of CSD between species and even within the same animal. Our model also explains 
differences between CSD in vivo and CSD in brain slices. 

Previous mathematical models of CSD have accounted for extracellular space and intracellular space 
(both somatic and glial). Here, we formulate a five-compartment continuum model for CSD (neurons 
consisting of two compartments, extracellular space consisting of two compartments, and vascular tree 
compartment), and perform numerical tests of the model under various conditions. A glial compartment 
is not considered in this model except in the context of potassium buffering. Finally, we discuss the 
results of simulations of the model and its relevant clinical and translational implications. 

2 Model 

To study the important new elements that affect and arc affected by CSD, we formulate a five-compartment 
continuum model. Neurons comprise two of the five compartments: a compartment representing the den- 
dritic processes (d), and a compartment representing the cell bodies (somatic compartment s). The 
extracellular space (e) also comprises two of the compartments, with the extracellular space proximal to 
cell bodies and proximal to dendrites each receiving a compartment. The vascular bed (y) comprises the 
remaining compartment. For computational tractability, we make the simplifying assumptions that were 
made in Yao et al. |16], that the intra- and extracellular volume fractions remain fixed. 

The neuronal and extracellular compartments include only the most-relevant ions (sodium, potassium, 
chloride) and channels that have been shown to be the most responsible for the instigation and spread of 
CSD [i"6Hl8| . In the somatic membranes, we include T-type sodium channels, P-type sodium channels, 
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delayed-rectifier potassium channels, A-type potassium channels, and the sodium-potassium ATPase. In 
the dendritic membranes, we additionally include NMDA channels. 

Our model assumes that the vascular compartment does not exchange fluid with the extracellular 
space, so that it does not affect the intra- and extracellular volume fractions. The effective diameters of 
proximal arterioles controls the blood flow rate. In turn, the vascular diameters are coupled to neuronal 
activity through extracellular potassium concentrations proximal to dendritic processes. 



2.1 Membrane potential and ion transport using a neuronal model 

The membrane potentials of the neuronal compartments £/ m> * (* is either s for somatic or d for dendritic), 
are governed by the coupled ordinary differential equations 

dendrite coupling 



Cm r; — — — Is.tot + nr > c2 {E m ,d — Em.s) (1) 

at £rLaO r] 



soma coupling 



Cm TT - — —Id,tat + r. n c2 {Em,s — E m ,d) (2) 

at ^tia"d 

where C m is the membrane capacitance per unit surface area (^farad/cm 2 ) for both the somatic and 
dendritic membranes, R a is the input resistance of the effective dendritic tree, 8d is the half length of the 
effective dendritic tree, and 7*,tot = J]] I*, ion, tot = J]] J]] I*, ion, channel are the spatially-dependent 

ions ions channels 

total cross-membrane ionic currents per unit surface area (mA/cm 2 ) for each neuronal compartment. 
Following Kager et al. [TTIHH] and Yao et al. [IB], the total cross-membrane currents, i*,tot ; are given for 
the three major ions (sodium, potassium, and chloride) and are the sum of the active and passive (leak) 
sodium and potassium currents, the chloride (leak) current, and the sodium-potassium exchange pump 
current (see Appendix |Aj . 

The local rates of change of the extracellular ions (Na + ,K + , and CP) are due to membrane ionic 
currents, diffusion of extracellular ions both within and between the extracellular compartments, and the 
buffering of extracellular potassium. They are given by 

through channels 



d([\on]e iS ) _ S s 

— r 7— it t ^s. ion, tot 

dt f e FV s 

diffusion within compartment 



JL( D . ^ 

~ I -is ion 



OX \ OX 

exchange with dendritic extracellular space 



Dion{Vd + Vs) „. , r. , , ,_x 
2^2y (l 10n k«J - [lOnJ e , s ) (3) 



d([ion] e , d ) _ S d T 

r\, r 7-1 t t d, ion. tot 

at feFVd 

d[ion} e , d 



^ I -^ion" 
OX \ OX 

Aon(Vs + Vd) n . , r. , v 

TjT^T ([ 10n \e, S ~ [lOn] M ) (4) 
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where S7*) are the surface areas of the neuronal compartments in the total fixed volume given by the 
sum of the fixed somatic V s , dendritic Vd, extracellular somatic f e V s , and extracellular dendritic f e Vd 
volumes. The equations for extracellular potassium are modified by adding the buffering flux given in 
Section 12.31 The equations for the rates of change of intracellular ions are 

through channels exchange with dendrites 



dt = -FV s Is ' io ^ ot + 2SlV s ([l ° nkd " [l ° llks) (5) 

d([ion] ij( j) S Aon(K + Vd) n . , [■• i \ 

j t = -py-Id,ion,tot + (l l0n ks - [ion]i, d ) . (6) 

Note that we have not included the diffusion of intracellular ions because any appreciable diffusion 
(relative to cell size) that would require the ions to first become an extracellular ion and then diffuse 
cxtracellularly. 

As CI - is the only extracellular anion in our simulations, to ensure electroneutrality, its initial extra- 
cellular concentrations ([Cl - ] ei *) were determined by the sum of the concentrations of the major cations: 
[Cl _ ] ej * = [Na + ] e> , + [K + ] e> ». We chose the initial intracellular concentrations of chloride ([Cl - ]^*) so 
that its Nernst potentials matched the membrane potentials. For electroneutrality intracellularly in the 
soma and dendrites, we added a non-specific immobile anion. Although calcium can have major effects on 
neurotransmitter dynamics and other secondary messenger effects, we ignored the calcium concentration 
because of its relatively small values. 



2.2 Exchange pumps 

During CSD, the ionic concentrations in the extra- and intracellular spaces are considerably displaced 
from steady-state. This displacement occurs primarily because of fluxes through voltage-gated sodium 
and potassium channels. Here, we include the sodium-potassium exchange pumps (Na-K-ATPase) in the 
neuronal membranes, whose primary role is to restore the ionic concentrations back to their homeostatic 
state. The ionic pumps are active and consume energy. When local oxygen levels are depleted, ATP is in 
short supply. Therefore, the function of the ionic pumps is related in our model to oxygen consumption 
and vascular flows. 

These pumps are involved in the movements of intracellular sodium and extracellular potassium 
against their electrochemical gradients require active ionic pumps that consume energy and the pumps 
are fueled by the dephosphorylation of ATP in the cell [TO] given by 

ATP + 3Na+ + 2K+ ADP + P t + 3Na+ + 2K+. 

When local oxygen in the tissue is very low, the normal ATP dynamics are perturbed. 

The Na + /K + -ATPase is a trans-membrane protein with two extracellular binding sites for potassium, 
three intracellular binding sites for sodium, and a single intracellular binding site for ATP. In each 
neuronal compartment, its potassium and sodium currents are given by pump = — 2/* jPU mp and 
4,Na,pum P = 3/* iPum p, respectively, where 

-^*,pump -^maxT*,pump,l'7*,pump,2 

and 

7*, P ump,i([ K+ k*, [Na + ]^„) 

= ( 1+ ^y 2 ( 1+ ^my 3 (7) 
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Figure 1. Oxygen availability affects the Na + /K + — ATPase. Shown is the relationship between 
tissue oxygen concentration and the oxygen-dependent portion of the pump rate (Eq[8]). Since some 
ATP is generated even in the absence of oxygen, the pump rate does not completely go to zero as the 
oxygen concentration approaches zero. The steady state oxygen concentration is 0.2mM. 



This expression is given in |17j and is a Hill equation for the pump rate under the assumption of plentiful 
ATP. We modify this equation so that it is valid even in the low- ATP regime by including another Hill-like 
term of the form (1 + [ATPjo/IATP])" 1 . The Hill-exponent is set to 1 to obey the known stoichiometric 
rule for the pump. 

To use this formulation for the pump activity, one needs to understand ATP dynamics within the cell. 
Cloutier et al. [20] provides a detailed system of ordinary differential equations for the generation of ATP 
through both aerobic and anaerobic respiration. One of their key assumptions is that adenylate kinase 
is in quasi-static equilibrium. This assumption is invalid for situations where the system is under high 
metabolic stress. Instead, we assume that production and consumption of ATP arc dominant, thereby 

finding a quasi-steady-state relationship between ATP and O2, ATP = consti + const2 ( ? 2 + o ) ■ 



Thus, we can describe the pump rate explicitly in terms of oxygen availability with the expression 



where a is the percentage (about 5%) of ATP production that is independent of oxygen [^T][22], the 
subscript denotes the equilibrium values, and [O2] is the tissue oxygen concentration (see Fig[T]). This 
expression indicates that the pumping rate will be reduced whenever there is a decrease of the oxygen 
level in the tissue. 

2.3 Potassium buffering 

Astrocytes play important roles in the instigation and propagation of CSD and in the neurovascular 
coupling through neurotransmitter- mediated signaling pathways [23]. A principal role of astrocytes is 
to buffer extracellular potassium [24] . This buffering is achieved through a variety of inward rectifying 
potassium channels in the glial membrane and is bolstered by the extreme polarity of glial cell membranes 
with membrane potential near the Nernst potential for potassium. For this study, we are not interested in 
the exact mechanisms of glial potassium buffering, but are interested in reproducing accurate potassium 
dynamics for our model. Thus, we incorporate astrocyte effects through empirical potassium buffering. 





(8) 
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Following Kager et al. [T7], we modeled the potassium-buffering flux, ^buffer, by the following first-order 
ordinary differential equation, valid at each position in the extracellular continuum 

"buffer = ~ = M + [K + ] e 3exp ( [K ^j; ~ 5 ' 5 ) - MA) - B), (9) 

where B (mM) is the free extracellular buffer concentration, the rate constants /i+ = /i_ = 8.0 x lCP 6 ms _1 
determine the speed at which potassium is buffered, and B = 200 mM is the effective total buffer con- 
centration. This equation describes strong buffering of extracellular potassium for concentrations above 
5.5 mM, but is limited by saturation of the finite buffer. As the amount of buffered potassium increases, 
re-release of potassium into the extracellular space becomes more favorable. The initial value of the free 
buffer concentration is set to maintain steady state when the extracellular potassium concentration is at 
its rest value (3.5 mM). 



2.4 Neurovascular coupling and oxygen supply 

We now describe how the tissue oxygen level is affected by and influences CSD. First, we assume that 
there exists an effective blood vessel radius r and that the cerebral blood flow rate (CBF) is given by 

r 4 

CBF = CBFo^r (10) 



r, 



in which the equilibrium values arc again denoted by the subscript 0. This expression is based on the 
empirical observation that blood flow through the small vessels where nutrient exchange primarily occurs 
can by modeled as Poiscuille flow where the volume flow rate is proportional to r 4 [25] . 

We use an empirical model for the effective vessel radius r, based on replicating the activity observed 
in many experimental studies on the subject |26U27j . Extracellular potassium is known to dilate vessels 
at lower concentration elevations (less than ~ 17 mM), and constrict vessels at higher concentration 
levels [28]. To reproduce this behavior, we assume that the effective vessel radius is given by 



constriction 



r([K + ] e ) = r exp 



[K+] e -3. 5 x 2 
a 

dilation 



(l + be -[([K+]e-10)/c]^ 

1 + 6 e -( 6 - 5 / c ) 2 



(11) 



which is a product of constricting and dilating terms. Since dendrites are long neuronal processes that 
occupy a greater volume than cell bodies, we take [K + ] e to be the concentration in the dendritic extracel- 
lular space. The parameter a (mM) controls the constriction response, b controls the amount of maximal 
dilation in the vessel, and c (mM) controls the dilation response. We found values of a = 50 mM, b = 0.18, 
c = 3mM to mimic the vascular response of Farr and David |28j (see Fig [2]) , who provide a model for 
the radius of cerebral arterioles based upon currents through potassium channels in the membranes of 
vascular cells. Using this expression gives us a simple relationship between the extracellular potassium 
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Figure 2. Relationship between vascular caliber and [K H 

function of extracellular potassium concentration [K + ] e . 



Effective vascular radius r as a 



concentration and CBF, namely, 

CBF = CBF 



exp 




1 + &e -[([K+] e -10mM)/c] : 
1 + 5g-(6.5mM/c) 2 



(12) 



We model the temporal evolution of the tissue oxygen concentration [O2], using a reaction-diffusion 
equation 



d[Q 2 

at 



= d 



d 2 [Q 2 
dx 2 



S 



with nonlinear source term given by 

[0 2 ] 6 - [0 2 ] 



S = CBF 



CBF 



[O a ] 6 - [O 2 ]o 

7pump,2 ([0 2 ]) 7pump,2 (0) 



7pump,2 ([0 2 ] 0) 7pump,2 (0) 



(1-7) 



+ 7 



7s, P ump,i([K + ] e>;i , [Na + ] z , s ) +7rf,pum P ,i([K + ] e|d , [Na + ]j, d ) 
27pum P a([K + ] e ,o,[Na + ] l , ) 



(13) 



The first term in the source is the amount of oxygen transferred from the blood stream to the tissue and 
is given by the product of CBF and the normalized concentration difference in oxygen tension between 
the blood and the tissue. The remaining terms represent the consumption of oxygen by the sodium- 
potassium exchange pump and other cellular processes that are assumed to remain steady during CSD. 
The pump consumption is given by the product of the equilibrium CBF and pump rate normalized by 
the steady-state pump rate. A fraction < 7 < 1 of the total oxygen consumption at steady-state is due 
to the pump. Experimental estimates of this parameter have ranged from as low as 0.10 [29] to as high 
as 0.70 [3D]. We include simulations over the full spectrum for thoroughness. Note that we have S = 
at steady-state. 



3 Numerical Tests 



In order to shorten numerical computational time, we sought to reduce the size of the computational 
domain to the minimal size required in order to achieve stable propagating waves that did not exhibit 
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boundary effects. In order to do this, we used a minimal amount of potassium to induce CSD. We 
found that for injections of Gaussian boluses with 120 micron width, a K + concentration of 15mM was 
sufficient to induce CSD. This concentration is near the experimentally-determined minimum threshold 
for CSD induction [31] . and confirmed to us that our phcnomcnological potassium buffer was behaving 
in a physiologically realistic manner. 

We performed one-dimensional simulations by breaking up our model into a system of ODEs solved 
using Matlab routine odel5s, with reflecting boundary conditions. The computational domain was 
set to a length of 2.7 cm, and discrctized into 46 grid points. This domain was sufficiently long for 
the propagating wave to become stable and for boundary effects to remain insignificant. For diffusion, 
we used a centered finite-difference approximation of the Laplacian operator, resulting in a system of 
ordinary differential equations for the ionic and oxygen concentrations at each grid point. We measured 
the duration of potassium elevation by taking the total amount of time that [K + ] Cit 2 is above 6 mM. The 
transient sodium channel had negligible effect on the ionic currents during CSD (data not shown), so we 
chose to omit it from these simulations. 

The variables in the dendritic and somatic compartments did not differ significantly in any of the sim- 
ulations. For this reason, we report results from the dendritic compartments and drop the corresponding 
prefix. 

3.1 Oxygen-clamped simulations 

First, we simulated our model with the oxygen coupling completely turned off, i.e., 7 = (see Eqll3j). 
Thus, the sodium-potassium- ATPase does not consume any oxygen, and we assume that the tissue 
oxygenation is held at its steady state value. In this situation, we obtain CSD waves with a velocity 
of 3.2 millimeters per minute and with extracellular potassium concentration increasing to a maximum 
value of 45.7 mM. The total duration of the elevated extracellular potassium concentration is 65.4 seconds 

(Fig E]). 

The potassium buffer acts similar in that in Kager ct al. [17j . It saturates rapidly, after which it is 
responsible for a net re-release of potassium into the extracellular environment. The buffer acts this way 
in all simulations as it is independent of the oxygen level in our model. For this reason, we omit further 
mention of the potassium buffer. 

3.2 Blood vessel-clamped simulations 

Next, we simulated CSD where the effective blood vessel diameter remained fixed, thereby fixing the 
maximal oxygen flux rate into the tissue. We varied the oxygen coupling parameter 7 between (no 
oxygen consumption by pump) and 1 (maximal oxygen coupling) in increments of 0.025. 

The results, shown in FigUJ illustrate the response of CSD to oxygen coupling. By increasing 7, we 
increase the duration and the amplitude of the CSD waveform. To get a deeper understanding of these 
results, Fig [5] shows the macroscopic observables. One sees that increases in 7 lead to greater drops 
in oxygenation. This deoxygenation results in decreases in the magnitude of the sodium-potassium- 
ATPase pump current, resulting in longer tails of recovery for both membrane potential and extracellular 
potassium concentration. 

3.3 Coupled vasculature 

From the previous results, one sees that oxygenation has a definite impact on susceptibility of the tissue 
to CSD and on its recovery. The next question is whether vasoconstriction has an impact on CSD. We 
performed a series of simulations where we activated the vascular coupling that we defined in Eqs. fTTI 
and 1131 We simulated an idealized vascular response, using parameters mimicking the vascular behavior 



9 




20 40 60 80 100 120 
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Figure 3. CSD in absence of oxygen consumption. Shown are the time-courses of the 
propagating waves of extracellular potassium concentration, free buffer capacity, membrane potential, 
and sodium-potassium- ATPase pump rate at positions 120 microns apart for simulations performed in 
the absence of oxygen consumption by the Na + /K + -ATPase. The Na + /K + -ATPase operates at rates 
determined solely by [K + ] e and [Na + ]i. 

of Farr and David [28]. In these simulations, we varied the oxygen coupling parameter, 7, in order to see 
the variability in behavior one might expect to see in an idealized CSD experiment. 

The results, shown in Fig [SJ illustrate the effects of vascular coupling. One can see that the effective 
blood vessel radius drops rapidly to approximately 40% of its original value. The drop leads to steep 
sustained reductions in the oxygen concentration during the metabolic challenge of the CSD wave. The 
deoxygenation is reflected in the pump rate, where a large increase of inward potassium current at the 
beginning of CSD is rapidly diminished as oxygen is depleted. The implications of this chain of events 
are visible in the extracellular potassium profile. 

As we noted earlier, vascular response is variable, depending on factors such species, individual, and 
metabolic states. In Fig [5J the effects of this variability are seen in simulations across a grid of values for 
the parameters a, 6, 7 (Eq. [TTj) . The duration of the CSD event decreases as the constriction parameter 
a is increased, and as the dilation parameter b is increased. 

In Fig[5J we have plotted the speed of CSD waves across the spectrum of possible vascular responses. 
These vascular responses are parameterized mainly by two parameters, the constriction parameter a, and 
the dilation parameter b (Eg. [TTj) . 
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Figure 4. Oxygen coupled to Na + /K + ATPase. Fixing the effective vascular radium while 
coupling oxygen to the sodium-potassium- ATPase gives a purely consumption-based view of oxygen 
during CSD. Defining coupling as the fraction of oxygen consumed by the pump at steady-state, these 
simulations show that all of the CSD characteristics - speed, duration, maximum [K + ] - increase with 
an increase in the coupling constant 7. The duration of the wave increases linearly from 7 = to 
7 « 0.4 before increasing linearly at a different rate when 7 > 0.4 

4 Discussion of major findings 

Since its discovery, CSD has been known to involve massive changes in vascular caliber, and hence, 
perfusion [6]. Up to this point, these changes, which can have profound effects on cortical function and 
thus on CSD itself, have not been incorporated into CSD models. Using our model, we were able to 
examine the effects of vascular activity on CSD. 

4.1 Relationship between vascular activity and CSD 

Integration of different vascular effects more closely approximates in vivo CSD and helps explain why CSD 
seems to be different in slice preparations. Our numerical experiments confirmed that oxygen delivery 
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Figure 5. Fixed vascular caliber. Time courses at a fixed position, 780/im downstream of the 
original stimulus. These simulations show the effects of oxygen consumption on CSD for a range of the 
oxygen coupling parameters, 7. Increasing 7 prolongs the duration and magnifies the amplitude of the 
CSD because it implies that the pump consumes more of the available oxygen, thereby resulting in 
larger oxygen depletions. 

plays a significant role in the dynamics of CSD. Even in the absence of constriction, oxygen depletion 
is seen. This fact, however, did not absolve the constriction from its implication in oxygen depletion. 
Additionally, vasoconstriction and vasodilation were seen to further modify the characteristics of the CSD 
wave, particularly in the low-7 regime. 

By varying the relationship between vasoconstriction and extracellular potassium, we achieved a con- 
tinuum of CSD responses that could explain real physiological variability across species, within members 
of the same species, and across repeated CSDs in the same animal. 

4.1.1 CSD uncoupled from vascular activity 

Our modeling of CSD uncoupled from oxygen concentration resembles current CSD models, which rely 
on similar assumptions and conductances. Such models are relevant to CSD in ex-vivo brain slices, 
where the perfusing medium has contact with an essentially unlimited source of oxygen. In this case, the 
pump is always able to perform at a rate determined by [K + ] e and [Na + ]i, and the total pump current is 
constrained only by the total number of available pumps present, provided that there is sufficient nutrient 
available. 
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Figure 6. "Typical" vascular coupling. Time courses at a fixed position, 780^m downstream of 
the original stimulus. Simulations performed using vessel of Farr and David |28j . After a short dilation 
period, the vessels constrict significantly to about 40% of their rest radii, before recovering. 



4.1.2 CSD coupled with vascular activity 

Coupling an inert vasculature into CSD, one obtains a sense of the metabolic demands that CSD places 
on the cortex. The large displacements of sodium and potassium require the use of energy to rectify, and 
the physical constraints placed on blood delivery by the physiologically reasonable assumption of finite 
volumc-flow-rate were seen to result in depletion of oxygen and reduction in potassium flux through the 
pumps. This finding is consistent with in vivo measurements of tissue oxygen and mitochondrial redox 
state even with mild dilation. Several studies have found that both variables moved into a more-reduced 
oxidative state [32]f33] . In our model, the depletion of oxygen is due to increased metabolic demand driven 
by the sodium-potassium-ATPase. 
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Figure 7. Effect of vascular coupling. When vasculature is coupled to [K ] e , the extra depletion of 
oxygen due to supply constraints causes the pump to operate at a slower rate than when 
vasoconstriction does not occur. As a consequence, the inward potassium flux through the pump is 
decreased. The result is prolonged recovery time. Shown are simulations where 7 = 0.7. 



4.1.3 The effects of vasoconstriction/vasodilation 

Furthermore, we show that vasoconstriction can both increase tissue's susceptibility to CSD and impair 
its ability to recover. These effects are due to the reductions in blood flow causing a significant additional 
drop in cortical oxygenation levels. In Fig [6j the vessel is seen to constrict to approximately 40% its 
original radius. This constriction is within the experimental range reported in Chang et al. |10j . Due to 
the power-law relationship between blood flow and effective blood vessel radius, even a small reduction 
in effective radius has large blood flow implications. A 60% constriction results in blood flow dropping to 
2.6% of its original value. This effect is visible in the oxygenation level, which undergoes further decreases. 
Our simulations show that both wave speed and recovery time increase when vessel constriction is greater. 

Our simulations also show that vasodilation plays a role in CSD. Predictably, vasodilation appears 
to precondition the tissue so that it is more able to withstand the increase of extracellular K + that 
accompanies the CSD wave. Due to the power law relationship, an increase of 14% in the vessel radius 
results in a 70% increase in blood volume flow rate. The result shown in Fig[9]is a tissue that is more able 
to withstand potassium elevations, thereby causing slow-downs in the CSD wave. Vasodilation (beyond 
the effective resting radius) seems to play a significant role in the recovery from CSD, causing decreases 
in the recovery time (Fig [5]) . 



14 




Figure 8. Duration of CSD for different vascular responses. Shown arc CSD durations plotted 
against percent dilation b (the maximum vasodilation), constriction parameter a (this parameter has 
the units of mM and is the width of a Gaussian curve that controls how fast r drops as [K + ] 
increases (Eq. fTTj) ), and oxygen coupling constant 7. Duration is defined by the length of time that 
potassium concentration is elevated to a level greater than 6 mM. Increasing a, (shown from 
30 — 80 mM) decreases the amount of constriction, resulting in quicker recovery from CSD. Likewise, 
increasing b, which increases the maximum dilation of the vessels, also reduces the duration of CSD. 




Figure 9. Speed of CSD for different vascular responses. Increasing the constriction parameter 
a makes the vessel less sensitive to potassium, weakening the resulting constriction. The speed 
decreases as a increases from 30 mM up to 80 mM, but increasing the amount of dilation by increasing 
b seems to have a greater effect. 



4.2 Clinical and translational implications 

If a major goal of CSD modeling is to understand its role in human disease, incorporation of perfusion 
and metabolism is an essential step. CSD is a near-maximal depolarization of brain tissue which of neces- 
sity challenges the ability of homcostatic mechanisms to compensate. CSD is also triggered by changes 
in perfusion or metabolism. Finally, during CSD, the relationships between brain activity, perfusion 
and metabolism - neurovascular coupling - are altered. These complex relationships can be expected 
to lead to a variety of responses, depending on the state of the brain tissue surrounding the depolar- 
ization. Biological data from animals and humans bears out this complexity and variability. Spreading 
depolarizations can be relatively innocuous - repetitive CSD in mouse over several days appears to cause 
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no overt injury [34) . However, they can also be quite harmful, enlarging infarct and contusion areas in 
both animals and humans [211253 . These deleterious effects are almost certainly due to alterations in the 
vascular response to tissue depolarization, and thus cannot be understood from a modeling standpoint 
without explicit incorporation of perfusion and metabolism. 

Our model more realistically represents conditions observed in the brain during experimental CSD 
and the spreading depolarizations of migraine and brain injury. Though it simplifies a complex vascu- 
lar/metabolic response, it has the distinct advantage of making specific, quantifiable predictions which 
can be used to generate hypotheses for further experimentation. A particular advantage is the ability 
to explore the whole "CSD/metabolic parameter space," which is not possible in biological experiments. 
This modeling could have important implications for the study of CSD in migraine, as the conditions 
which could generate such a massive depolarization in an awake behaving person remain obscure. 

4.3 Assumptions, limitations, and future directions 

In this study, in order to make our model as widely applicable as possible, we have not considered the 
geometry of any particular vascular network. Our lumped model is assumed to be a good approximation 
of oxygen delivery dynamics in the brain, in general, for gray-matter. Further insight into the system 
may be gained by targeting specific tissue types for simulation. 

Our simplified model relating arterial effective diameter to K + is likely applicable during CSD but does 
not mechanistically explain all the subtleties of neurovascular coupling under more normal conditions. 
Though K + is involved in the coupling of neural and astrocytic activity to blood flow [36] , other mediators 
including arachidonic acid derivatives, purines, nitric oxide, and possibly neurotransmitters are involved 
as well [7JH31I3Z]- For the purposes of CSD, however, with its supra-physiological swings in extracellular 
potassium, K + likely provides a good leading-order approximation of vascular behavior. Many of these 
other mediators also work via their effects on potassium channels |28j . 

Swelling of individual cells [17l|38l|39] and the tissue as a whole [10] appears to occur during CSD. 
Whether cell-swelling plays an important role in CSD is a controversial topic, as recent computational 
studies have cast some doubt on its impact |TB]13D]. For this reason, we have omitted cell swelling from 
our model. 

Finally, this model is not equipped to account for the complexity of the post-CSD state. Chang et 
al. [10] found that extracellular K + levels remained constant during a long-lasting hypoperfusion and 
depolarization that was found to follow CSD. As this phase may be clinically relevant, future modeling 
will focus on trying to understand the K + -independent mechanisms involved in the etiology of this period. 
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A Cross-Membrane Currents and Parameter Values 

We assume that the membrane currents in the soma used by Kager et al. [T7] apply here. The total cross- 
membrane current in the soma is given by the sum of the active and passive sodium, potassium, chloride, 
and nonspecific ionic currents. The sodium current is / s ,Na,tot = / s ,Na,p + ^s,Na.Lcak + ^s,Na,Pump where 
In&.p is the persistent sodium current, I Si Na,ioak is the sodium leak current, and I s ,Na,pump is the current 
through the pump. The active and passive potassium currents are I Sj K,tot = 4,k,dr + Is,k,a + Is,K,\eak + 
^s,K.pump where I s .k,dr is the potassium delayed rectifier current, I s .k.a is the transient potassium 
current, is,K,pump is the potassium current through the pump, and / Si K,ieak is the potassium leak current. 
The passive chloride leak current is is, leak- 

In the dendritic compartment, the NMDA currents, id.Na,NMDA and id,K,NMDA, are added to the total 
current. Note that we have removed I s ,Na,T, the fast transient sodium current, since this current was 
shown in |16j to not make any fundamental difference in the way CSD propagates. The mathematical 
expressions for these channels do not differ across the dendritic and somatic compartments so we omit 
the compartment prefix (s,d) from this point forward. 

The cross-membrane currents are modeled using: i) the Goldman-Hodgkin-Katz (GHK) formulas for 
the active membrane currents given by 



ion. GHK 



m p h q - 



9ion,GHKFE m 


[ion] i - cxp 




[ion] e 





l-exp(-^ 







(14) 

which is suitable when there is a large difference in concentrations between the ICS and ECS compart- 
ments, as argued by Koch and Segev and ii) the Hodgkin-Huxley (HH) model for the passive leak 
currents given by 

hon,HH = 9ion,HH (E m — £ion) ■ (15) 

In these general expressions for the GHK and HH types of currents, <7ion,GHK and <7ion,HH are the con- 
ductances associated with the channels for ion = Na + and K + , and m and h are the activation and 
inactivation gating variables, respectively, for the different GHK-modelcd channels that are ion-specific. 
Note that the (7i on ,HH conductances for the sodium and potassium leak currents and the <?i ea k conductance 
for the chloride leak current are assumed to be constant. The variables, E- lon , are the Nernst potentials 
for ion = Na + and K + given by 

£ ion = 01ogM£ (16) 
[ionJ, 

For the chloride leak current, the equivalent Nernst potential is — 70mV. The factors in the parameter 
cf> = RT/F are R, the universal gas constant, T = 310K, the absolute temperature, and F, the Faraday 
constant. 

The gating variables, m and h, satisfy the following relaxation equations 

^ = a m (l -m) - /3 m m, (17) 
at 

— = a h (l-h)-p h h, (18) 

and the values of a and f3 are given in Table [TJ along with the exponents p and q [T5] . The extracellular 
volume is assumed to be 15% of the intracellular volume, i.e., V e — 0.15V;. Other relevant parameter 
values are: R = 8.31 (mV coulomb/mM K) and F = 96.485 coulomb/mM. 
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We use the following procedure to choose the initial (equilibrium) values for the gating variables and 
ion concentrations. We first set the membrane potential at E. m = —70 mV and the sodium and potassium 
concentrations (listed in Table [5]) , from which we compute the parameters a and /3 (using the formulas 
in Table [lj . We then compute the equilibrium values of m and h as 

a m au 

Tfl = ft = . 

a m + /3 m ' ah + fih 

Next we choose the leak conductances, <?Na,ieak and <?K,icakj by setting 

^Na,tot = ^K.tot = 0. 

Finally, we determine the chloride leak conductance in Ji ea k by assuming that gi ca k = lO^Na.ieak- 
Table 1. Parameter values for active membrane ionic currents, from [P7 ] [T8 ] . 



Currents 
mA / cm 2 


9ion,GHK 

S/cm 2 


Gates 
m p h q 


Voltage-Dependent Rate Constants 




2 xl0~ 6 


m 2 h 


l 

m 6(l+cxp[-(0.143E m mV- 1 +5.67)]) 
o exp[-(0.143E m mV- 1 +5.67)] 
^ m 6(l+exp[-(0.143£; m mV- 1 +5.67)]) 

a h = 5.12 x 10- 8 exp[-(0.056 J E m mV~ 1 + 2.94)] 

a . i.6xi0" 6 

>^ h l+exp[-(0.2B m mV- 1 +8] 


^K,DR 


10 xlO" 5 


o 

Tfl 


m m u - ulu i„ cxp [_(o.2£; m mV- 1 +6.98)] 

/3 m = 0.25exp[-(0.25 J E m mV- 1 + 1.25)] 


^K,A 


1 xlO" 5 


m 2 h 


n . _nno E m mV- 1 +56.9 

«m u.uz, 1 _ exp [_( . 1 _ BmmV -i +5<69 )] 

a _ n m 7^ S m mV- 1 +29.9 

U - U - L ' U exp(0.1E m mV- 1 +2.99)-l 

a h = 0.016 expl-^.OSe^mV" 1 + 4.61)] 

a . 0.5 

l Jh l+expI-tO^BmmV- 1 +11.98)] 


^NMDA 


1 xlO" 5 


mh 


am - i+ex P ( i3 - s -[ 1 + 4 ir M - i ) 

f3 m = 0.5 - a m 

^ ~ 2000(l+exp[f K+ ' em M 7 ; i - 6 - 75 ]) 

ft = 5x 10" 5 - a h 



The initial resting ion concentrations consistent with those in [TT1IT8] are obtained by modifying 
the sodium-potassium exchange pump function and running the model equations until a steady state is 
reached. These values are given in Table [2] together with other parameter values. 



Table 2. Initial resting values and other relevant parameter values for the computations, following 
Kager et al. QUEE]. 



Parameter 


Value 


Unit 


R a (input resistance of dendritic tree) 


1.83 x 10 b 


ohms 


d s (diameter of soma) 


5.45 x 1(T 4 


cm 


dd (diameter of dendrite) 


0.939 x 10~ 4 


cm 


5d (half-length of dendrite) 


4.5 x 1CT 2 


cm 


A s (surface area of soma) 


1586 x 1(T S 


cm 2 


Ad (surface area of dendrite) 


26732 x 10~ 8 


cm 2 


V s (volume of soma) 


2160 x 10" 12 


cm 3 


Vd (volume of dendrite) 


5614 x 10~ 12 


cm 3 


C'm (membrane capacitance) 


7.5 x lO" 5 


s /fl cm 2 


I mo .x (Na+/K+-ATPasc rate) 


1.48 x 10~ 3 


mA / cm 2 


E m 


-70 


mV 


[K + ] e 


3.5 


mM 




133.5 


mM 


[Na+] e 


140 


mM 


[Na+]. t 


10 


mM 


[O 2 ] 


2 x 10~ 2 


mM 


CBF 


2.5 x 10~ 2 


mM/s 


Do, 


0.5 x 10~ 5 


cm 2 /s 


[0 2 ] 6 


4 x 10~ 2 


mM 



